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We show that the classical Rosenbluth method for sampling self-avoiding walks [2, 0] can be 
extended to a general algorithm for sampling many families of objects, including self-avoiding poly- 
gons. The implementation relies on an elementary move which is a generalisation of kinetic growth; 
rather than only appending edges to the endpoint, edges may be inserted at any vertex providing 
the resultingobjects still lie within the same family. We implement this method using pruning and 
enrichment [3| to sample self-avoiding walks and polygons. The algorithm can be further extended 
by mixing it with length-preserving moves, such pivots and crank-shaft moves. 
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Monte Carlo simulations of self-avoiding walks (SAWs) 
and self-avoiding polygons (SAPs) on regular lattices are 
a major tool for the study of polymer statistics Q ■ While 
kinetic growth algorithms [J HI have been used to sample 
SAWs to great success, it is unclear how they might be 
applied to SAPs. SAPs are models of ring polymers, 
plasmids (mitochondrial DNA) [f| 0, 0] and appear in 
the zero-component limit of the N- vector model @. In 
this paper, we generalise a growth algorithm for SAWs 
and show how it may be used to sample SAPs and other 
objects. 

The Rosenbluth method for sampling self-avoiding 
walks (SAWs) is a classical algorithm dating back to the 
1950s P, 0. This method found new application with 
the development of a pruned and enriched implementa- 
tion, called PERM, due to Grassberger Q. This was 
in turn further extended using flat-histogram and mul- 
ticanonical methods 0, [l(| • These new algorithms have 
found man y a pplications in the modelling of polymers 

(eg.PElll). 

The Rosenbluth method samples SAWs by growing 
conformations from a single vertex and iteratively adding 
edges from the final vertex to its unoccupied nearest 
neighbours. Let u> be a SAW of n edges and let rrik be 
the number of unoccupied nearest neighbour vertices of 
its endpoint when it is truncated after k edges (after the 
k th iteration). The probability of sampling u> is given by 



SAWs of length n is 



PrH = n 
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If mfe = the conformation is trapped and a new con- 
formation is started. This attrition makes it difficult to 
sample long SAWs and sampling by this method is not 
uniform. The sample bias can be removed by weighting 
each SAW by W{u) = Pr(w) -1 . The mean weight of all 



W(u)Pi(w) 



(2) 



where c n is the number of SAWs of length n. Thus one 
can estimate the number of conformations by computing 
the average weight of sampled SAWs. Sample averages of 
observables are computed by taking weighted averages. 

The original Rosenbluth algorithm converges very 
slowly because of attrition and the large variance in 
weights of longer SAWs. Grassberger's pruned and en- 
riched implementation Q constituted a major advance 
that allows the algorithm to reduce these effects. This 
significantly improves its efficiency and applicability. 

In this paper we demonstrate how the Rosenbluth al- 
gorithm can be generalised; SAWs can be sampled by 
inserting edges at any vertex rather than only at the end- 
point. This idea may be extended to a growth algorithm 
for SAPs and other lattice objects. We also show how to 
combine it with length preserving moves such as pivots. 

Atmospheres 

Let uj be a SAW on the square lattice, starting from 
the origin. A positive endpoint atmospheric edge of u 
is an edge on the lattice that can be appended to the 
last vertex of u> to extend its length by one while re- 
specting self- avoidance (see Figured]). The size of the 
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A SAW and its two positive endpoint atmospheric 
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positive atmosphere of to, a + (w), is the number of posi- 
tive atmospheric edges. We abuse our notation and use 
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a + to denote both the set of positive atmospheric edges 
and its size. Adding positive endpoint atmospheric edges 
increases the length of the SAW. 

Appending a positive atmospheric edge to ui to ob- 
tain uj' creates a linkage (see Figure [5]). Deleting 
the last edge from uj' gives uj; we define the negative end- 
point atmosphere of uj' to be this edge. We denote the 
size of the negative atmosphere of uj' by a_ (a/) and we 
again abuse notation by also using this symbol to denote 
the set of such edges. In the present context a_(a/) = 1, 
but below we consider more general positive and negative 
atmospheres. 
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FIG. 2: A schematic picture of the linkages between Cn 
and Cn+i- There is one conformation with empty positive 
atmosphere. 

By counting the number of linkages, we see that 

#linkages = a + (ui) = a_(u/) = c n+1 . (3) 



This implies that 
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where a_(o/) = 1 and the averages are taken over the 
uniform distribution. 

This observation can be used to estimate growth 
constants and free-energies of SAWs and bond trees 
EH, EH • We extend these definitions and show how 
they lead to a significant generalisation of the Rosenbluth 
algorithm. 

Define the positive atmosphere, a+(w), to be the num- 
ber of ways that an edge can be inserted into uj at any 
of its vertices so that a SAW is obtained (see Figure [3]). 
Note that there are SAWs with empty positive atmo- 
sphere. 
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FIG. 3: The SAW on the right is obtained from the SAW on 
the left by inserting a north edge at the black vertex. This 
is one of its eleven positive atmospheric edges. It has three 
negative atmospheric edges. 



Inserting a positive atmospheric edge into u> results in 
a new SAW, uj'. The negative atmosphere, a_(o/), is the 
number of ways that an edge can be deleted to obtain 
a SAW. Linkages are created as above and equation 
holds by the same arguments. Simulations show that 
the distribution of atmospheres are narrowly peaked (see 
Figure HJ . 

Generalised atmospheric Rosenbluth method 

The generalised atmospheres, a±, can be used to define 
a generalised Rosenbluth method for sampling SAWs. 
The algorithm starts with a single vertex, <pg , and grows 
a sequence of SAWs, tp = tpo, . . . , tp n , by inserting a pos- 
itive atmospheric edge at each iteration. The confor- 
mation ifk+i is obtained from tpk by inserting an edge 
chosen uniformly from the available positive atmospheric 
edges, a+((fk)- We call this the Generalised Atmospheric 
Rosenbluth Method (GARM). This method generalises 
the percolation based algorithms for trees in [171 EH • 

A sequence of n + 1 SAWs, tp = ipo, • ■ ■ , tp n , is obtained 
after n iterations with probability 



Pi(tp\tp ) = JJa+(yfe-i)" 



(5) 
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Since a given conformation can be obtained in several 
different ways, this is not the probability of obtaining 
the last SAW in the sequence. As such we give a weight 
to the sequence of SAWS, not only to compensate for 
the non-uniform sampling probability, but also to take 
into account this degeneracy. The weight of a sequence 
of SAWs, tp = po, . . . , p n is 



w(tp) = ii 



-{pk-i 
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a-(pk) 



if n > 1 and W(tp) = 1 if n = 0. 

The mean weight of sequences of length n + 1 is 

(W) n =Y,W{tp)?r{tp\tp ) = c n . 



(6) 



(7) 



To see this, consider all the sequences that end in a par- 
ticular SAW t of length n. It suffices to show that 



i = j2 w &) p r(^ i = j2 n a -^y 
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where the sums are over all sequences that end in t. Note 
that one must consider all the possible choices of atmo- 
spheric edges, so that there are n\ sequences that end in 
the SAW made up of n east edges. 

We reinterpret the product of negative atmospheres as 
the probability of returning to the single vertex under 
the following process. Starting at tp n , we delete negative 
atmospheric edges from pk to obtain pk-i iteratively. 
The probability of realising po along the sequence tp is 



PrM^) = Y[Pr(p, 
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Since all sequences that end at r must return to <po by this 
process, summing over ip gives equation ([5]) as required. 
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FIG. 4: Typical evolution of W{ip) 1/n for both Rosenbluth 
and GARM sampling. Shown are 20 samples from each. 
While GARM does suffer from attrition, SAWs sampled by 
Rosenbluth have larger variance and a higher rate of attri- 
tion. The insets show the distributions of positive and nega- 
tive atmospheres per vertex (top) and endpoint atmospheres 
(bottom); the peak heights have been normalised to 1. 

The above proof becomes trivial in the case of the end- 
point atmosphere since a_ = 1 and each SAW is obtained 
in exactly one way. The proof breaks down in models in 
which a given conformation cannot be reached by insert- 
ing positive atmospheric edges. 

In Figure [5] we show that data obtained by a pruned 
and enriched implementation of GARM for SAWs agrees 
with exact enumeration data from [19(. We generated 
SAWs of a maximum of 200 edges with approximately 
10 6 trajectories consisting of a total of 2.5 x 10 8 samples. 
This took about an hour on a laptop computer. This 
algorithm for the square-lattice generalises to SAWs on 
any graph with finite maximal degree. 

Extensions to polygons, trees and animals 

We now extend this algorithm to SAPs on the square 
lattice. In a previous paper, we defined the positive at- 
mosphere to be the locations in which a single edge can 
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FIG. 5: A plot of estimated from GARM data up to 

length 200. Exact enumeration data up to length 71 taken 
from [l9d is shown for comparison. 



be replaced by a n conformation of three edges and the 
negative atmosphere was defined by the inverse of this 
process [2(| ■ This definition is insufficient for GARM 
since there are many conformations that are not obtain- 
able from the unit square; for example the 2x2 square. 

We generalise the notion of positive atmospheres of 
SAPs by considering all the pairs of vertices at which 
anti-parallel edges may be inserted to obtain a longer 
SAP. The negative atmosphere is defined by finding all 
pairs of edges that may be removed to obtain a SAP. See 
Figure [SI All SAPs have non-zero positive atmosphere. 
Since the atmospheres now consist of pairs of edges we 
have 

{a+)2n _ P2n+2 
(a-)2n+2 P2n 

where pi n is the number of SAPs of length 2n. 



FIG. 6: A SAP and the insertion of a pair of anti-parallel 
edges. This SAP has positive atmosphere of 21 and negative 
atmosphere of 4. 

In Figure [7] we show that data obtained by a PERM- 
like implementation of GARM for SAPs agrees with ex- 
act enumeration data from [2l[ . We generated SAPs of a 
maximum of 200 edges with approximately 4 x 10 5 tra- 
jectories consisting of a total of 5 x 10 7 samples. This 
took a few hours on a laptop computer. 

This algorithm does not simply generalise to three di- 
mensions with this definition of atmospheres. If ipa is 
chosen to be a unit square, then knotted conformations 
cannot be reached since inserting atmospheric edges does 
not allow strand passages. 

GARM can be applied to lattice bond trees. The algo- 
rithm is then closely related to the algorithm in [Tt], [3 • 
In this case we define the positive atmosphere by looking 
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FIG. 7: A plot of pl/ n as estimated from GARM data up to 
length 200 for even n. Exact enumeration data up to length 
110 taken from [2ll ] is shown for comparison. 
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at all the vertices at which an edge can be inserted to 
obtain a valid tree. When inserting an edge at a given 
vertex, one must be careful to consider all the possible 
ways of distributing the incident branches between both 
ends of the new edge. One may similarly define posi- 
tive and negative atmospheres for animals. The positive 
atmosphere is defined by all the ways in which an edge 
may be inserted at vertices; unlike the tree case, some in- 
serted edges will create cycles and so are double counted 
as they can be inserted from either vertex. The negative 
atmosphere is defined by the inverse of this process and 
atmospheric edges that are not cut-edges will be double 
counted. This algorithm works for site-trees; this may 
be easily implemented by defining positive and negative 
atmospheres in terms of leaves, but more general atmo- 
spheres are possible. 

The implementation of the GARM algorithm requires 
rapid calculation of the positive and negative atmo- 
spheres. For SAWs the positive and negative atmo- 
spheres are 0(n) while for SAPs they are 0(n 2 ). At 
present we are able to compute the atmospheres in 0(n) 
time for SAWs and 0(n 2 ) for SAPs. Since the atmo- 
spheres must be computed at each iteration, the time to 
produce a conformation of length n is 0(n 2 ) and 0(n 3 ) 
for SAWs and SAPs respectively. 

Conclusions 

The definitions of atmospheres above were limited to 
positive and negative since they either increase or de- 
crease the number of edges. We can generalise this 
further by including the notion of neutral atmospheric 
moves, ao, which change the conformation without 
changing its size — for example a pivot move. At each 
iteration the algorithm chooses uniformly to add an edge 



from the positive atmosphere or to apply a neutral atmo- 
spheric move. The probability of obtaining a sequence tp 



M-i 

Pr(tp\ipo) = Yl {a+(fk-i) + aoiifk-i))^ 1 ■ (11) 

k=l 

and the corresponding weight is 



w{ V ) = n 

k=l 



q + ((/?fc_i) +a (y>fc-i) 
a-(tp k ) + ao(<pk) 
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where \tp\ is the number of conformations in the sequence 
tp. The average weight of all sequences ending in a confor- 
mation of size n is c n ; the proof is as above. This addition 
makes it possible to sample SAPs in three dimensions and 
higher since the pivot algorithm is ergodic [221 ]. 

The algorithm can also be adapted to include Boltz- 
mann factors (as per 0, [23J) so as to compute free en- 
ergies. Further extensions such as multicanonical or flat 
histogram methods, such as those developed in [§, |T(| 
are possible. We are currently investigating techniques 
to compute atmospheres more efficiently as this will im- 
prove the convergence of the GARM algorithm. 
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